nonlin_multi_eqn_mult_var.f90

Source Code

module nonlin_multi_eqn_mult_var
    use iso_fortran_env
    use nonlin_types
    use nonlin_error_handling
    implicit none
    private
    public :: vecfcn
    public :: jacobianfcn
    public :: vecfcn_helper
    public :: equation_solver
    public :: nonlin_solver

    interface
        subroutine vecfcn(x, f, args)
            !! Describes an M-element vector-valued function of N-variables.
            use, intrinsic :: iso_fortran_env, only : real64
            real(real64), intent(in), dimension(:) :: x
                !! An N-element array containing the independent variables.
            real(real64), intent(out), dimension(:) :: f
                !! An M-element array that, on output, contains the values
                !! of the M functions.
            class(*), intent(inout), optional :: args
                !! An optional argument to allow the user to communicate with
                !! the routine.
        end subroutine

        subroutine jacobianfcn(x, jac, args)
            !! Describes a routine capable of computing the Jacobian matrix
            !! of M functions of N unknowns.
            use, intrinsic :: iso_fortran_env, only : real64
            real(real64), intent(in), dimension(:) :: x
                !! An N-element array containing the independent variables.
            real(real64), intent(out), dimension(:,:) :: jac
                !! An M-by-N matrix where the Jacobian will be written.
            class(*), intent(inout), optional :: args
                !! An optional argument to allow the user to communicate with
                !! the routine.
        end subroutine
    end interface

    type vecfcn_helper
        !! Defines a type capable of encapsulating a system of nonlinear
        !! equations of the form: F(X) = 0.  This type is used to establish the
        !! system of equations to solve, and provides a means for computing
        !! the Jacobian matrix for the system of equations, and any other
        !! ancillary operations that may be needed by the solver.
        procedure(vecfcn), private, pointer, nopass :: m_fcn => null()
            !! A pointer to the target vecfcn routine.
        procedure(jacobianfcn), private, pointer, nopass :: m_jac => null()
            !! A pointer to the jacobian routine - null if no routine is 
            !! supplied.
        integer(int32), private :: m_nfcn = 0
            !! The number of functions in m_fcn.
        integer(int32), private :: m_nvar = 0
            !! The number of variables in m_fcn.
    contains
        procedure, public :: set_fcn => vfh_set_fcn
        procedure, public :: set_jacobian => vfh_set_jac
        procedure, public :: is_fcn_defined => vfh_is_fcn_defined
        procedure, public :: is_jacobian_defined => vfh_is_jac_defined
        procedure, public :: fcn => vfh_fcn
        procedure, public :: jacobian => vfh_jac_fcn
        procedure, public :: get_equation_count => vfh_get_nfcn
        procedure, public :: get_variable_count => vfh_get_nvar
    end type

    type, abstract :: equation_solver
        !! A base class for various solvers of nonlinear systems of equations.
        integer(int32), private :: m_maxEval = 100
            !! The maximum number of function evaluations allowed per solve.
        real(real64), private :: m_fcnTol = 1.0d-8
            !! The convergence criteria on function values.
        real(real64), private :: m_xtol = 1.0d-12
            !! The convergence criteria on change in variable values.
        real(real64), private :: m_gtol = 1.0d-12
            !! The convergence criteria for the slope of the gradient vector.
        logical, private :: m_printStatus = .false.
            !! Set to true to print iteration status; else, false.
    contains
        procedure, public :: get_max_fcn_evals => es_get_max_eval
        procedure, public :: set_max_fcn_evals => es_set_max_eval
        procedure, public :: get_fcn_tolerance => es_get_fcn_tol
        procedure, public :: set_fcn_tolerance => es_set_fcn_tol
        procedure, public :: get_var_tolerance => es_get_var_tol
        procedure, public :: set_var_tolerance => es_set_var_tol
        procedure, public :: get_gradient_tolerance => es_get_grad_tol
        procedure, public :: set_gradient_tolerance => es_set_grad_tol
        procedure, public :: get_print_status => es_get_print_status
        procedure, public :: set_print_status => es_set_print_status
        procedure(nonlin_solver), deferred, public, pass :: solve
    end type

    interface
        subroutine nonlin_solver(this, fcn, x, fvec, ib, args)
            !! Describes the interface of a nonlinear equation solver.
            use, intrinsic :: iso_fortran_env, only : real64
            use nonlin_types, only : iteration_behavior
            import equation_solver
            import vecfcn_helper
            class(equation_solver), intent(inout) :: this
                !! The [[equation_solver]]-based object.
            class(vecfcn_helper), intent(in) :: fcn
                !! The [[vecfcn_helper]] object containing the equations to
                !! solve.
            real(real64), intent(inout), dimension(:) :: x
                !! On input, an N-element array containing an initial estimate 
                !! to the solution.  On output, the updated solution estimate.
                !! N is the number of variables.
            real(real64), intent(out), dimension(:) :: fvec
                !! An M-element array that, on output, will contain the values 
                !! of each equation as evaluated at the variable values given 
                !! in x.
            type(iteration_behavior), optional :: ib
                !! An optional output, that if provided, allows the caller to 
                !! obtain iteration performance statistics.
            class(*), intent(inout), optional :: args
                !! An optional argument to allow the user to communicate with
                !! the routine.
        end subroutine
    end interface

contains
! ******************************************************************************
! VECFCN_HELPER
! ------------------------------------------------------------------------------
    subroutine vfh_set_fcn(this, fcn, nfcn, nvar)
        !! Establishes a pointer to the routine containing the system of
        !! equations to solve.
        class(vecfcn_helper), intent(inout) :: this
            !! The [[vecfcn_helper]] object.
        procedure(vecfcn), intent(in), pointer :: fcn
            !! The function pointer.
        integer(int32), intent(in) :: nfcn
            !! The number of functions.
        integer(int32), intent(in) :: nvar
            !! The number of variables.
        this%m_fcn => fcn
        this%m_nfcn = nfcn
        this%m_nvar = nvar
    end subroutine

! ------------------------------------------------------------------------------
    subroutine vfh_set_jac(this, jac)
        !! Establishes a pointer to the routine for computing the
        !! Jacobian matrix of the system of equations.  If no routine is
        !! defined, the Jacobian matrix will be computed numerically (this is
        !! the default state).
        class(vecfcn_helper), intent(inout) :: this
            !! The [[vecfcn_helper]] object.
        procedure(jacobianfcn), intent(in), pointer :: jac
            !! The function pointer.
        this%m_jac => jac
    end subroutine

! ------------------------------------------------------------------------------
    function vfh_is_fcn_defined(this) result(x)
        !! Tests if the pointer to the subroutine containing the system
        !! of equations to solve has been assigned.
        class(vecfcn_helper), intent(in) :: this
            !! The [[vecfcn_helper]] object.
        logical :: x
            !! Returns true if the pointer has been assigned; else, false.
        x = associated(this%m_fcn)
    end function

! ------------------------------------------------------------------------------
    function vfh_is_jac_defined(this) result(x)
        !! Tests if the pointer to the Jacobian calculation routine has been
        !! defined.
        class(vecfcn_helper), intent(in) :: this
            !! The [[vecfcn_helper]] object.
        logical :: x
            !! Returns true if the pointer has been assigned; else, false.
        x = associated(this%m_jac)
    end function

! ------------------------------------------------------------------------------
    subroutine vfh_fcn(this, x, f, args)
        !! Executes the routine containing the system of equations to
        !! solve.  No action is taken if the pointer to the subroutine has not
        !! been defined.
        class(vecfcn_helper), intent(in) :: this
            !! The [[vecfcn_helper]] object.
        real(real64), intent(in), dimension(:) :: x
            !! An N-element array containing the independent variables.
        real(real64), intent(out), dimension(:) :: f
            !! An M-element array that, on output, contains the values
            !! of the M functions.
        class(*), intent(inout), optional :: args
                !! An optional argument to allow the user to communicate with
                !! the routine.
        if (this%is_fcn_defined()) then
            call this%m_fcn(x, f, args)
        end if
    end subroutine

! ------------------------------------------------------------------------------
    subroutine vfh_jac_fcn(this, x, jac, fv, args)
        !! Executes the routine containing the Jacobian matrix if
        !! supplied.  If not supplied, the Jacobian is computed via finite
        !! differences.
        class(vecfcn_helper), intent(in) :: this
            !! The [[vecfcn_helper]] object.
        real(real64), intent(inout), dimension(:) :: x
            !! An N-element array containing the independent variables defining 
            !! the point about which the derivatives will be calculated.
        real(real64), intent(out), dimension(:,:) :: jac
            !! An M-by-N matrix where, on output, the Jacobian will
            !! be written.
        real(real64), intent(in), dimension(:), optional, target :: fv
            !! An optional M-element array containing the function values at x. 
            !! If not supplied, the function values are computed at x.
        class(*), intent(inout), optional :: args
            !! An optional argument to allow the user to communicate with
            !! the routine.

        ! Parameters
        real(real64), parameter :: zero = 0.0d0

        ! Local Variables
        integer(int32) :: j, m, n, lwork, flag
        real(real64) :: eps, epsmch, h, temp
        real(real64), pointer, dimension(:) :: fptr, f1ptr
        real(real64), allocatable, target, dimension(:) :: wrk

        ! Initialization
        m = this%get_equation_count()
        n = this%get_variable_count()

        ! Input Checking
        flag = 0
        if (size(x) /= n) then
            flag = 2
        else if (size(jac, 1) /= m .or. size(jac, 2) /= n) then
            flag = 3
        end if
        if (flag /= 0) error stop flag

        ! Process
        if (.not.this%is_fcn_defined()) error stop NL_UNDEFINED_FUNCTION_ERROR
        if (associated(this%m_jac)) then
            ! Call the user-defined Jacobian routine
            call this%m_jac(x, jac, args)
        else
            ! Compute the Jacobian via finite differences
            if (present(fv)) then
                lwork = m
            else
                lwork = 2 * m
            end if

            ! Local Memory Allocation
            allocate(wrk(lwork))
            f1ptr => wrk(1:m)
            if (present(fv)) then
                fptr => fv(1:m)
            else
                fptr => wrk(m+1:2*m)
                call this%fcn(x, fptr, args)
            end if

            ! Establish step size factors
            epsmch = epsilon(epsmch)
            eps = sqrt(epsmch)

            ! Compute the derivatives via finite differences
            do j = 1, n
                temp = x(j)
                h = eps * abs(temp)
                if (h == zero) h = eps
                x(j) = temp + h
                call this%fcn(x, f1ptr, args)
                x(j) = temp
                jac(:,j) = (f1ptr - fptr) / h
            end do
        end if
    end subroutine

! ------------------------------------------------------------------------------
    function vfh_get_nfcn(this) result(n)
        !! Gets the number of equations in this system.
        class(vecfcn_helper), intent(in) :: this
            !! The [[vecfcn_helper]] object.
        integer(int32) :: n
            !! The function count.
        n = this%m_nfcn
    end function

! ------------------------------------------------------------------------------
    function vfh_get_nvar(this) result(n)
        !! Gets the number of variables in this system.
        class(vecfcn_helper), intent(in) :: this
            !! The [[vecfcn_helper]] object.
        integer(int32) :: n
            !! The number of variables.
        n = this%m_nvar
    end function

! ******************************************************************************
! EQUATION_SOLVER
! ------------------------------------------------------------------------------
    pure function es_get_max_eval(this) result(n)
        !! Gets the maximum number of function evaluations allowed during
        !! a single solve.
        class(equation_solver), intent(in) :: this
            !! The [[equation_solver]] object.
        integer(int32) :: n
            !! The maximum number of function evaluations.
        n = this%m_maxEval
    end function

! --------------------
    subroutine es_set_max_eval(this, n)
        !! Sets the maximum number of function evaluations allowed during
        !! a single solve.
        class(equation_solver), intent(inout) :: this
            !! The [[equation_solver]] object.
        integer(int32), intent(in) :: n
            !! The maximum number of function evaluations.
        this%m_maxEval = n
    end subroutine

! ------------------------------------------------------------------------------
    pure function es_get_fcn_tol(this) result(x)
        !! Gets the convergence on function value tolerance.
        class(equation_solver), intent(in) :: this
            !! The [[equation_solver]] object.
        real(real64) :: x
            !! The tolerance value.
        x = this%m_fcnTol
    end function

! --------------------
    subroutine es_set_fcn_tol(this, x)
        !! Sets the convergence on function value tolerance.
        class(equation_solver), intent(inout) :: this
            !! The [[equation_solver]] object.
        real(real64), intent(in) :: x
            !! The tolerance value.
        this%m_fcnTol = x
    end subroutine

! ------------------------------------------------------------------------------
    pure function es_get_var_tol(this) result(x)
        !! Gets the convergence on change in variable tolerance.
        class(equation_solver), intent(in) :: this
            !! The [[equation_solver]] object.
        real(real64) :: x
            !! The tolerance value.
        x = this%m_xtol
    end function

! --------------------
    subroutine es_set_var_tol(this, x)
        !! Sets the convergence on change in variable tolerance.
        class(equation_solver), intent(inout) :: this
            !! The [[equation_solver]] object.
        real(real64), intent(in) :: x
            !! The tolerance value.
        this%m_xtol = x
    end subroutine

! ------------------------------------------------------------------------------
    pure function es_get_grad_tol(this) result(x)
        !! Gets the convergence on slope of the gradient vector
        !! tolerance.
        class(equation_solver), intent(in) :: this
            !! The [[equation_solver]] object.
        real(real64) :: x
            !! The tolerance value.
        x = this%m_gtol
    end function

! --------------------
    subroutine es_set_grad_tol(this, x)
        !! Sets the convergence on slope of the gradient vector tolerance.
        class(equation_solver), intent(inout) :: this
            !! The [[equation_solver]] object.
        real(real64), intent(in) :: x
            !! The tolerance value.
        this%m_gtol = x
    end subroutine

! ------------------------------------------------------------------------------
    pure function es_get_print_status(this) result(x)
        !! Gets a logical value determining if iteration status should be
        !! printed.
        class(equation_solver), intent(in) :: this
            !! The [[equation_solver]] object.
        logical :: x
            !! True if the iteration status should be printed; else, false.
        x = this%m_printStatus
    end function

! --------------------
    subroutine es_set_print_status(this, x)
        !! Sets a logical value determining if iteration status should be
        !! printed.
        class(equation_solver), intent(inout) :: this
            !! The [[equation_solver]] object.
        logical, intent(in) :: x
            !! True if the iteration status should be printed; else, false.
        this%m_printStatus = x
    end subroutine

! ------------------------------------------------------------------------------
end module
nonlin_multi_eqn_mult_var.f90 Source File

Contents

Modules

Source Code

Source Code
